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Abstract 

The transition from hquid metal to silicate rock in the cores of the terrestrial 
planets is likely to be accompanied by a gradient in the composition of the outer 
core liquid. The electrical conductivity of a volatile enriched liquid alloy can be 
substantially lower than a light-element-depleted fluid found close to the inner 
core boundary. In this paper, we investigate the effect of radially variable elec- 
trical conductivity on planetary dynamo action using an electrical conductivity 
that decreases exponentially as a function of radius. We find that numerical 
solutions with continuous, radially outward decreasing electrical conductivity 
profiles result in strongly modified flow and magnetic field dynamics, compared 
to solutions with homogeneous electrical conductivity. The force balances at the 
top of the simulated fluid determine the overall character of the flow. The rela- 
tionship between Coriolis and Lorentz forces near the outer boundary controls 
the flow and magnetic field intensity and morphology of the system. Our results 
imply that a low conductivity layer near the top of Mercury's liquid outer core 
is consistent with its weak magnetic field. 

Key words: Variable electrical conductivity, numerical dynamos, Geodynamo, 
Mercury 



1. Introduction 

Variations in the physical properties of fluids in planetary dynamos define 
the character of the observed intrinsic magnetic field (e.g., strength, geometry 
and time variability). Changes in the electrical conductivity of the fluid as a 
function of depth may become relevant in the context of terrestrial and gas giant 
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planets. In this paper we explore (with a focus on the terrestrial planets) how 
the radial variation of electrical conductivity in planetary cores may result in 
changes to dynamo-generated magnetic fields. 

1.1. Terrestrial planets 

The cores of the terrestrial planets are composed principally of iron, with 
minor but significant amounts of nickel and lighter elements. It has long been 
known that an iron-nickel core would have too high a density to be compati- 



ble with Earth's moment of inertia and seismic data (e.g., Birch, 1952 Poirier 



1994). A compatible Earth core density model can result from the inclusion of 
about 8% by weight of one or more light elements. Detailed models of core com- 
position are based primarily on the constraints of seismology, mineral physics, 
geochemistry, metallurgy, and cosmochemisty. Silicon, sulphur and oxygen are 
the primary candidates for the light elements. Sulphur is likely to be a signifi- 
cant component but the depletion of light elements in the process of accretion 
in the inner solar system limits sulphur to about 2 wt% in the core. 

Recent reviews of core differentiation and composition distinguish between 
models considering Silicon versus those considering Oxygen as the primary light 
element. Composition models give weight percents of Fe 85%-88%, Ni 5% 
Si 0-7%, O - 0-4% and S 2% (e.g .McDonough[ [20031 [Wood et aH[2QQ6| . 
Solidification of a more or less pure iron-nickel inner core may exclude the lighter 
elements, which would then be enriched in the outer core. The ratio of the inner 
core radius (1221 km) to the core radius (3480 km) is 0.35 and the mass of the 
inner core is only about 5% of the total core mass. So the bulk composition of 
the outer core is only slightly different than that of the whole core. 

Convection in the liquid outer core is driven by a combination of compo- 
sitional and/or thermal buoyancy. Thermal buoyancy available to drive con- 
vection and dynamo action originates primarily from the latent heat of solid- 
ification at the inner core boundary (ICB) and possibly from cooling at the 
core-mantle boundary (CMB). Secular cooling at the CMB does not guarantee 
a source of convective instability since heat can be conducted through a stably 
stratified layer (a super-adiabatic heat fiux is needed). Compositional buoyancy 
originates at the ICB due to light element enrichment in the residual liquid as- 
sociated with inner core solidification. A source of compositional buoyancy near 
the CMB could come from precipitation from a silicate enriched layer at the top 
of the core, leaving a light element depleted, heavy residual liquid and a silicate 



sediment layer at the CMB (Buffett et al. 2000). 



Temperature and pressure dependence of electrical conductivity may lead 



to sizable variations with depth. (Stacey and Anderson, 2001). Assuming a 



well- mixed outer core and a Si concentration of X^. = 0.25 in a dual species 
alloy, the authors found the most extreme variation in the fiuid core of Earth 
results in a factor of 1/2 difference for the electrical conductivity and about 3/4 
for the thermal conductivity from the ICB to the CMB. Although the effects of 
different impurities in the alloy have not been explored experimentally, [Stacey 
and Anderson (2001) argue that the effects of impurities other than Si do not 



significantly change their estimates. A later revision (Stacey and Loper 2007) 
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found that the variation in electrical and thermal conductivities had been over- 
estimated due to the assumption of treating Fe as an electronically simple metal. 
Their later results predict a difference in the electrical conductivity of a factor 
of 0.78 between that of the CMB when compared to the ICB's. Those estimates 
were made under the assumption of a well- mixed outer core. However gradients 
in material properties would be amplified in a compositionally stratified layer. 

The existence of a thermally and/or compositionally stably stratified layer 
near the top of the Earth's core has been both, suggested and argued (iJacobs 



1975 Fearn and Loper, 1981 Lister and Buffett 1998 Braginsky, 1993| 2007). 



Here we will focus on the possibility of a compositionally stratified layer near the 
top of the core since the electrical conductivity in such a layer would decrease 
with radius (due to the higher light-element concentration when compared to 
the bottom of the fiuid core). Two mechanisms by which a compositionally 
stratified layer can grow are by chemical diffusion from the mantle directly 



to the top of the core (Lister and Buffett, 1998) and by buoyant transport 



of light element enriched residual liquid from inner core solidification ( MoffattI 



and Loper 1994). Chemical diffusion is a slow process and would not result 
in a layer thickness greater than about 10 km (Lister and Buffettl 1998). On 



the other hand, buoyant rise of light element rich residual liquid could be an 
efficient process to build a layer of greater thickness. A thick layer with a 
significantly anomalous density would likely be detected as a seismically fast 
layer. Seismological constraints have been so far inconclusive. However, recent 
seismological results using outermost core waveforms seem to be consistent with 
the existence of a low density layer of about 100 km thickness ( [Alexandrakis 



and Eaton, 2007 Tanaka, 2007) 



It has been confirmed that Mercury has a liquid iron core (Margot et al. 



2007). Furthermore, it is likely that Mercury's weak intrinsic magnetic field is 



generated by a dynamo in its liquid outer core. Neither the details of Mercury's 



core composition, nor the size of its inner core is known ( |Hauck et al. , 2004 



Solomon et al. , 2008 Heimpel and Kabin, 2008). Observations of Mercury's 



contractional lobate crustal faults imply a planetary radial contraction of about 
3 km, small compared to an estimated 17 km of contraction that would result 
from a completely solidified core. This implies that Mercury may have a rel- 
atively large outer core, perhaps earth-like proportion. For a thick-shell outer 



core to exist in Mercury, given its small size, Hauck et al. ( 2004, ) estimates that 
the light elements (such as S or Si) in the core have a relatively high concen- 
tration. This suggests the possibility of a large and stratified Hermean outer 
core with a thick, and stably stratified outermost layer. Such a layer would be 
compositionally stratified (in contrast to thermally stratified models previously 



proposed, e.g., Christensen[ 2006) with a radial increase in the proportion of a 
light elements, and a radial decrease in electrical conductivity. 



1.2. Gas giant planets 

The magnetic field in the gas giants is generated within a metallic hydrogen 
region. Experimental results have found that the transition from metallic to 
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molecular hydrogen yields a wide range of pressures where the electrical con- 
ductivity is non-negligible thus varying slowly with depth in planetary interiors 
(Nellis, 2000). The internal structure has been deduced to first order based 



on measurements of the moment of inertia and total mass of each planet (e.g., 



Guillot 2005). However, constraints on hydrogen and helium mixtures at high 



pressures need to be found in order to better determine the internal structure 
of the gas giants. The additional effect of helium may complicate further the 
underdetermined internal layering of gas giants (e.g., Stevenson, 2008). The 



metallic to molecular hydrogen transition is of particular interest from the point 
of view of the internal dynamics. The depth and radial extent of this transition 
are important for magnetic field generation and in understanding the observable 
magnetic field morphology. This is a task for future investigation. 

In this paper, we present results from a numerical dynamo model with ra- 



dially varying electrical conductivity (Gomez-Perez, 2007). To implement the 



radially variable conductivity, we modified an existent numerical code (origi- 
nally Magic 2.0, |Wicht[ |2QQ2| , that uses the Boussinesq approximation. We 
focus this study on the effect of the varying conductivity on the dynamo action, 
and on the generated magnetic field. With this new implementation we per- 
formed a set of tests to analyze its consistency with previously published work. 
In section 12] we include a review of numerical models that worked with strati- 
fied liquids. In section [3] we present the necessary modification to the dynamo 
equations, and to the numerical code. The parameters explored for twenty runs 
studied in this paper are included in section |4j We present the results in sec- 
tion [5) the discussion and conclusions are found in section |6] and [?) respectively. 
We also included a table of symbols in appendix [Al 



2. Stratified planetary interiors in numerical simulations 

The study of planetary and solar dynamos has evolved rapidly in the past two 
decades due to code development and the accessibility to powerful computers 
used to solve numerically the equations of motion of rotating fluids. Solutions 
of numerical dynamos share characteristics with Earth's dynamo with respect 



to temporal variability, evolution and mean geometry (e.g., Kageyama et al 
Glatzmaier and Roberts[ |1995[ [Kuang and Bloxham[ |1997| [Christensen 



1993^ 



et al. 



1998). More recently, numerical models have been successful in reproduc- 



ing non-dipolar and weak fields such as those observed for the ice giants (e.g., 
Gomez-Perez and Heimpel 2007 Stanley and Bloxham, 2006). Nevertheless, 



these models of planetary dynamos rely on strong assumptions that do not nec- 
essarily represent realistic physical conditions expected for planetary interiors. 
For example, due to hardware limitations, numerical models cannot simultane- 
ously resolve planetary- and small-scale flow, both of which are prevalent for low 
viscosity fluids such as liquid iron or metallic hydrogen in a planetary setting. 

2.1. Buoyancy stratification 

Christensen and Wicht (2008), and Christensen (2006[) implemented a stable 



stratified (non-convecting) layer at the top of the electrical conductive fluid in 
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planetary dynamo models. They showed that small scale magnetic fields, may 
be generated in the deep interior of the electrically conductive fluid, shielded 
by a stably stratified layer. This outer layer acts as a filter that mainly damps 
the rapidly varying small scale components by the magnetic skin effect. They 
conclude that the axisymmetric field of Saturn and the relatively weak field of 
Mercury may be explained by a stably stratified layer at the top of the dynamo 
region. Similarly, Stanley and Mohammadi (2008) using buoyantly stratified 



fluids, found that zonal winds developed in a stable layer can affect the magnetic 
field time evolution rendering a secular variation of thermally stratified models 
inconsistent with Earth's observed field. 



2.2. Density stratification 

We discuss briefly the effects of density variation (anelastic versus Boussi- 
nesq approximations) in this section, since they are closely linked to those of 
the electrical conductivity. In order to solve the anelastic dynamo equations 



numerically, Glatzmaier (1984) used a poloidal and toroidal decomposition for 



the velocity and magnetic fields (for the decomposition of solenoidal vectors 
the reader may refer to Chandrasekhar 1961) . This allows for the use of a 



pseudo-spectral algorithm, which has been extensively used in planetary and so- 
lar dynamo simulations (e.g., Glatzmaierl 1984[|Dormy et aL , 1998[ |Ch ristensen| 



et al. , 1999; Sakuraba and Kono^ 1999; Clune et al. , 1999). This mathematical 



approach requires field vectors to be solenoidal (i.e., divergence free) and thus, 
approximations are used for the continuity equation, \/ - {pu) -\- d p / dt = 0, where 
u is the velocity, p the density of the fluid and t is time. If the time scale of con- 
vection is large compared to the time for sound waves to travel across the depth 
of the shell, then dp/dt <C V- (pu) and the anelastic approximation, V- (pu) = 0, 
is obtained. In addition, if the density changes in the fluid are much smaller 
than the averaged density, one obtains the Boussinesq approximation, V-u = 0, 
where the buoyancy is determined by temperature changes exclusively, but the 
density perturbations with respect to the mean are neglected. 

Differences between the results obtained using Boussinesq and anelastic ap- 
proximations have been compared in the context of magnetic convection in 2D 
models (Evonuk and Glatzmaierl 2004). They find significant effects due to 



strong density stratification (about three orders of magnitude variation in den- 
sity). The length scales of the top and bottom thermal plumes differ due to 
the compressibility of the liquid. 3D effects on the flow are very significant as 
well, and although the effect on the density variation is likely to be important in 
three-dimensional systems, a detailed analysis of density variations in 3D flows 
needs to be completed. 



Anufriev et al. (2005) compared Boussinesq and anelastic approximations 



analytically and found that they differ in the thermodynamic equations. In the 
latter, the rate of work to expand or contract the fluid is included in the energy 
balance, while the Boussinesq treats a thermal balance as the total energy bal- 
ance. The anelastic approximation takes into account a more complete picture 
by including cooling or heating done by the flow. They propose a change in the 
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thermal boundary conditions in Boussinesq models in order to allow for a realis- 
tic energy balance in the system. This suggestion should be seriously considered 
in future numerical simulations that use the Boussinesq approximation. 

We are not aware of previous numerical simulations that studied the effects 
of varying electrical conductivities in the context of planetary interiors. Anelas- 
tic models of the interior of stars include variations on magnetic and viscous 
diffusivities, and thermal expansion coefficient changing as a function of back- 
ground density, ex p~^l^ ( Feat herst one et al. 2007 Clune et al. 1999). Their 



numerical models include high resolution runs, specially targeted for stellar dy- 
namos (e.g., the heat transport throughout the simulated volume is set up to 
follow hydrogen burning at the deep interior). Here, we study the effects of 
radially varying electrical conductivity on the flow and on the magnetic field 
generation in the liquid cores of terrestrial planets. 

3. Methodology 

S.l. Governing equations 

To describe the behavior of an electrically conductive fluid bounded by co- 
rotating spherical shells, we use the magnetohydrodynamic equations in the 
frame of the rotating fluid. Consider a spherical shell with inner and outer 
boundary radii and respectively, rotating with an angular velocity ft. 
Using the Boussinesq approximation, these equations are written as: 

^ + (u • V)u - V^u j + 2z X u 

= -VP+^^fT+-i-(VxB) xB, (1) 

Pr Qo Pm"" 

V-u = 0, V-B = 0, (2) 

r^T 1 

^+u.VT = -V^T, (3) 

dB „ , 1 

'dt 

where u and B are the velocity and magnetic induction vectors, respectively; t 
is time; T and P are the temperature and pressure scalars, respectively; g is the 
acceleration of gravity which changes linearly with depth, and go is the value 
of this acceleration at the outer boundary; r is the radial unit vector, and z is 
the unit vector in the direction of the angular momentum of the rotating frame. 
^(^) — is normalized magnetic diffusivity, where A(r) and = X{ri) 
are the magnetic diffusivities, in units of m^s~^, of the liquid core as a function 
of r and of the solid inner core respectively. 

Equations [l] to |4] are expressed in terms of the following non-dimensional 
parameters: the Rayleigh number. 



= V X (u X B) - -pi^V X (a(V X B)) , (4) 



Ra=^^^^^ (5) 
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where a is the thermal expansion coefficient and hi is the thermal diffusivity. 
The Ekman number, 



which is the ratio between viscous and Coriolis forces in the system. The Prandtl 
number, 



Pr=-, (7) 

which is the ratio between the viscous and the thermal diffusivities. And a 
modified magnetic Prandtl number. 



Pm* = — , 

A 7 



(8) 



which is the ratio between the viscous diffusivity of the fluid and the magnetic 
diffusivity at the inner boundary, A(r^) = A^. Note that the magnetic Prandtl 
number changes as a function of radius, we used the maximum magnetic Prandtl 
number for the definition of the non-dimensional units used in our numerical 
model. 

The second term in the right hand side of equation |4] may be written as: 



Pm* 



V X 



(a(V X B)) 



1 



Pm* 



-AV^B + VA X (V X B) 



(9) 



in which the variable conductivity modifies the term VA x (V x B) exclusively. 

The magnetic induction equation is modified by the spatial variability of 
A(r). Since B is solenoidal, we define the magnetic field as a function of two 
scalar potentials (see Chandrasekhar , 1961 Appendix III): 



B = V X V X + V X , 



(10) 



where $ is the toroidal and ^ the poloidal potential. From this decomposition 



the time variation of the potentials (Gomez-Perez, 2007): 



= f-(Vx(uxB)) 



A L 



H 



f • V X (V X (u X B)) 



A L 



H 



Lh 



Pm* 



dr ' 



(11) 



(12) 
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where Vr = r • V. It is convenient to write the horizontal component of the 
angular momentum operator as Lh = —v'^Vh^ with V h = V — Vr- The term 
— AV^B in equation [9] introduces a zero order effect of a varying electrical dif- 
fusivity, it results on a decrease of the ohmic dissipation with depth, for both, 
poloidal and toroidal components of the magnetic field. The radial derivative 
of the diffusivity has a first order effect, involving exclusively the poloidal com- 
ponent. This is controlled by the term VA x (V x B) in equation [9] 

3.2. Radially variable electrical conductivity 

The radial change in electrical conductivity, a{r) = , where fiQ is 

the magnetic permeability of vacuum, is modeled with a piecewise continuous 
function whose derivative is also continuous. The outermost part decreases 
exponentially while the interior changes as a power of the radius. One may 
write the normalized electrical conductivity as: 




^{r-rm) 



(13) 



where is the inner boundary, is the radius where were the electrical con- 
ductivity changes from a polynomial to an exponential function, a determines 
how fast the exponential function decreases outside the conductive volume, am 
is chosen to be the value of the normalized conductivity where the functions 
match, a{rm) = cfm^ and k — a{am — l)(^m — (see figure [l]). We choose 

to implement a continuous differentiable function in order to have an analytical 
solution for A and VA. This function for the electrical conductivity is chosen to 
model a layering between an light-element-enriched fluid with a light-element- 
poor fluid. During the solidification of the core, heavy elements are preferen- 
tially deposited on the inner core. As a result, light element enriched liquid 
rises buoyantly to the top of the outer core, just below the CMB. A more real- 
istic representation of the electrical conductivity variability, and of miscibility 
of iron and light element enriched iron at high pressures, needs to be studied. 
From an experimental point of view such measurements are very challenging. 
First principles numerical simulations (e.g., "Caracas and Verstraet^ 2QQ9| ) may 



lead to some useful results and a more realistic characterization of A(r) may be 
implemented in the future. 



4. Dynamo regimes studied 

We study self-sustained dynamos in two main regimes. First, we include two 
sets of cases with relatively low i?a, where the resultant fields are non-reversing 
and dipole dominated. Second, we analyzed strongly forced and time variable 
cases where dynamos result in dipolar magnetic fields with strong multipolar 
components. We chose a values for Ra = S.ARac = 2.5 x 10^ (set A), Ra = 
7.8Rac = 8.5 x 10^ (set B) and Ra = 16.8Rac = 1.3 x 10^ (sets C and C), 
where the value of Rac was calculated using the empirical formula for the onset 
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Planetary radius, r/r 



Figure 1: Normalized electrical conductivity as a function of radius. The electrical conduc- 
tivity is normalized with the conductivity of the inner core, rm was chosen to be a grid point 
and thus Xm = 0.91 ^ 0.9 and Xm = 0.82 ^ 0.8 are labelled in the text as Xm = 0.9 and 
Xm = 0.8 respectively. 



of convection in rotating spherical shells from Al-Shamali et al. (2004). For sets 
A, C and C, E = 10~^ and E = 10~^ for set B. Common parameters used for all 
sets are a Prandtl number Pr = 1, and a radius ratio x — ~ =0.35. For sets A, 
B and C non-slip and fixed temperature boundary conditions are used at both, 
bottom and top boundaries. We also include cases with top and bottom free- 
slip boundaries, i.e., set C. For this set, all other parameters identical to those 
of set C. The models account for an electrically conductive inner core where 
electrical conductivity is homogeneous in the inner core, continuous through 
the ICB, and may vary radially in the outer core. The electrical conductivity 
for the mantle was chosen to be zero. All simulations in sets A, C and C have 
the same resolution, with 61 radial levels for the outer core and 17 for the inner 
core; 288 {Imax = 120) grid points in the azimuthal direction and half of that in 
the latitudinal direction. For set B the latitudinal and azimuthal resolution was 
increased to 256 and 512 respectively {Imax = 170) and 65 radial levels were 
used. Hyperdyffusivities were not used for any of the runs presented in this 
paper. For the purpose of normalization only, the magnetic Prandtl number, 
Pm* = 5, is defined by the ratio between the viscous diffusivity of the fiuid, z^, 
and electrical diffusivity of the inner core, (which is the minimum diffusivity 
in the simulated volume). 

For each set A, B, C and C we compare two cases with homogeneous elec- 
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trical conductivity in the fluid with cases with electrical conductivity varying 
as a function of radius. To label these cases Xm = ^m/i^o is used, where is 



the same used in equation [13) Cases with homogeneous conductivity are: non- 
magnetic case (from here on referred to as Xm = 0) and a case with constant 
electrical conductivity (xm = oo). The radially variable electrical conductivity 
cases A, B, C and C have Xm ^ 0.8, Xm ^ 0-9 and Xm ^ 0.98, see flgure [l] 
The shape of the electrical conductivity function is deflned by equation [13] For 
all cases a = 10, and am = h 



5. Results 

Signiflcant changes in the internal dynamics of the simulations are obtained 
by varying the electrical conductivity proflle. Energy time series for sets A and 
C are included in flgures|2]and[3) The time variability of the solution varies with 
Xm- We flnd quasi-stationary solutions for set A for Xm = 0.8 and Xm = 0.9. 
The solutions for sets B, C and C did not result in quasi-stationary solutions, 
as may be seen for set C in flgure |3] 

We include Table [T] with time-averaged values, where the average has been 
taken over two viscous diffusion times. The standard deviation from the time- 
averaged value is also included. 

The Reynolds number. Re = ^^"^^ , for Xm = is always greater than for the 
corresponding Xm = oo. In particular, the toroidal kinetic energy is decreased 
due to the presence of Lorentz forces. 

Values of the averaged Elsasser number, A = ^ . ^ in the fluid outer core 
and at the top of the simulated fluid are included in columns four and flve in 
Table [l] It is important to note that the deflnition of A uses instead of A(r) 
to be consistent with the normalization and the non-dimensional parameters. 
Using the total magnetic energy, M, the axisymmetric dipolar magnetic energy, 
M^^ and the total and axisymmetric energy at the CMB, M^, and M^, respec- 
tively; we calculate the ratio of the axisymmetric dipole to the whole magnetic 
energy for the whole fluid (Q^ = M^/M) and at the CMB (Q^ = M^/Mo). 
These values are shown in columns 7 and 6, respectively. 

The axisymmetric dipole component is dominant for all runs in set A, see 
in column six in Table [l] In contrast, the magnetic fleld energy is distributed 
over higher multipoles in sets B, C, and C, due to the higher Ra/Rac value. 
For these three sets, the variable conductivity cases exhibit broader length scale 
features and weaker flelds at the CMB than case Xm = co, this is due to the 
diffusion of magnetic fleld through the low electrical conductivity volume. 

Observable characteristics of dynamos are limited to the radial component 
of the magnetic fleld at the CMB and its temporal variability. It is important to 
understand how the flow dynamics is correlated to the surface fleld. Quantities 
such as in Table [l] serve as an indication of the observable magnetic fleld. 
In contrast, characterizes how the magnetic fleld being carried by the core 
flow behaves. We compare how much of the fluid core axisymmetry is carried 
to the CMB surface fleld. To do this, we take the ratio of the normalized 
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Figure 2: For set A, time series of the total kinetic energy in the fluid (in back) and the total 
magnetic energy (in grey). From top to bottom increasing values of Xm] Xm = 0, Xm = 0.8, 
Xm = 0.9, Xm = 0.98, and Xm = oo in rows one to flve. 



surface axisymmetric dipole energy and the normalized fluid core axisymmetric 
dipole energy, this quantity is shown as a function of Xm in figure |4] With 
the exception of set A a tendency for decreasing Q^/Q^ is found for increasing 
Xm- This is reasonable since higher harmonics are preferentially damped by the 
diffusion through the non-conducting volume. Thus, models with larger non- 
conductive volumes result in stronger axisymmetric dipoles at the top surface. 
The exception to this in set A is given by the average of Xm = 0.9, the quasi- 



12 



10^ 
10^ 
104 

10^ 
10^ 
104 

10^ 
10^ 
104 

10^ 







0.5 



1.5 







0.5 



1.5 



u « i M>!'! i" !«ii j j i i f >!'i it * i < " ! o <*> :~ ! i i :iii*« i * ' ' i' MXjMu >« n^»''*^: i( ii n i . i . i i^:'* ' :ti i L->y j f I " im. 'v ,~.".'O ii j riij 







0.5 







0-) 


10^ 




W 






104 




10^ 












10^ 






w 






104 



1.5 







0.5 



1 



1.5 







0.5 1 

Viscous time 



1.5 



Figure 3: Similar to figure |2] but for set C. x-n 
Xm = oo in rows one to five. 



3, Xm = 0.8, Xm = 0.9, Xn 



0.98, and 



stationary solution results in an increased axisymmetry when compared to Xm = 
0.8. It is possible that this is a consequence of the periodic variation (see panel 
3 in figure [2| reinforcing and stabilizing the axisymmetric components of the 
magnetic field. 

In figure|5]equatorial profiles of the temperature field are shown. A behaviour 
change is evident between the homogeneous non-magnetic and magnetic solu- 
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Figure 4: Time average of the ratio of the relative axisymmetric dipole energy at the CMB 
(Qo) and the relative axisymmetric dipole energy over the fluid core (Q^). The ratio has 
been calculated for all runs where magnetic flelds are non-zero and the color refers to sets A 
(Black), B (dark grey), and C (light grey). Values of simulations in set C are shown with red 
squares. 

tions, with the variable conductivity solutions as intermediate states. For set 
A, columns one to three (xm = 0, O-S, 0.9) exhibit a similar behaviour. The 
presence of a mean zonal flow bending the hot thermal plumes in the prograde 
direction is visible, while for the fifth column {xm = oo) there is no evidence 
for a strong axisymmetric zonal flow, and the plumes have a dominant radial 
direction. Column four (xm = 0.98) shows a weak prograde tilt. For set B, a 
transition is seen between Xm = 0.9 and Xm = 0.98, although there is also a sig- 
nificant morphology difference between the non-magnetic case and the variable 
conductivity x = 0.8 and Xm = 0.9. For the third row (set C) there is a differ- 
ence between the finite Xm runs and Xm = oo in the characteristic azimuthal 
wave number of the plumes. In Xm = oo of set C there are six distinguishable 
hot plumes while for all others the wave number is as high as nine. 

In figure |6j snapshots (at the same instants as for figure |5| of equatorial 
profiles of the axial vorticity are shown. Corresponding equatorial profiles of 
|B| in figure [t] are included for all sets. 

In set A there is a clear correlation between the cyclonic vortices (in black) 
with the magnetic field magnitude maxima (in white). Similarly to the tempera- 
ture, the wave number of the vorticity profiles of the variable conductivity cases 
is equal to that of the non-magnetic case and allows for smaller scale features. 
For set A, as well as seen in the temperature profiles, the presence of zonal flow 
bends the vortices in a prograde direction. Runs with Xm = 0.8 and Xm = 0.9 
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mimic this prograde tendency and the magnetic field is regularly organized in 
azimuth with a dominant wave number of six. Runs Xm = and Xm = 0.98, 
do not show a clear prograde tilt with increasing radius, and have a higher dom- 
inant wave number than runs with lower values of Xm- The correlation between 
cyclonic vortices and magnetic field maxima can be inferred for sets B, C and 
C, but it is not as evident as for set A. The length scale of the magnetic field 
decreases with decreasing Xm- 

The radial component of the magnetic field at the top of the simulated fluid 
is shown in figure [S] 

The symmetry of the magnetic field refiects the symmetry of the internal 
fiow. Similarly, as found by |Heimpel et al.| ( |2QQ5| ), for quasi-geostrophic fiows 
(such as those in set A) the thermal plumes generate Taylor columns in the fiow 
that transport magnetic field lines, resulting in a magnetic field morphology that 
correlates well with the flow and the thermal structure. For the weakly forced 
system (set A), the case with a homogeneous electrical conductivity results 
in a dipolar dominated dynamo with low non-axisymmetric components. The 
variable conductivity cases, where the magnetic field generated by the strong 
dynamo region has diffused through the low electrical conductivity fluid, exhibit 
a dominant dipolar fleld as well, which is highly axisymmetric (more so than the 
homogeneous case). In contrast, the strongly forced systems (sets B, C and C), 
with a more disorganized flow regime, result in non-axisymmetric flelds with 
signiflcant higher degree components. 

There is an identiflable transition in flow and magnetic fleld morphology that 
corresponds to different values of the conductivity proflle radius ratio Xm- For 
cases with thick variable conductivity layers (low Xm), thermal and flow flelds 
resemble those of non-magnetic convection (see flgures 5 and 6). In contrast, for 
relatively thin variable conductivity layers (high Xm) result in flow and magnetic 
flelds that closely resemble the homogenous conductivity (xm = oo) cases. Such 
transition may be identifled in Table [l] as a large increase in the CMB Elsasser 
number, A^^^, for increasing Xm- 

This transition in magnetic fleld morphology and magnitude corresponds as 
well to a change in the relative strength of Coriolis and Lorentz forces. Time 
averages of the Coriolis and Lorentz forces as a function of radius for sets B 
and C are shown in flgures |9] and 10 The volume averages has been taken over 
cylinders of radius r that have a symmetry axis parallel to z. In set B, runs 
where A^^^ is large, the top of the core exhibits a balance between Lorentz and 
Coriolis forces, with the exception of the region included in the Ekman layer 
(i.e., the last three grid points at the top of the simulated volume). For set 
C, the Coriolis is balanced by Lorentz in Xm = For both sets of magnetic 
cases with low A^^^, the sum of Lorentz, Coriolis and Reynolds stresses is 
completely dominated by the Coriolis force, which is balanced by viscous forces 
at the top of the core. In this set for the case with Xm = 0.98, it is hard to asses 
whether there is a balance between Coriolis and Lorentz forces or whether the 
Coriolis dominates. 
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Figure 9: Time average of magnetic (blue) and Coriolis (red) forces as well as Reynolds stresses 
(green) over concentric cylinders at radii from to Tq. In panels top to bottom increasing 
values of Xm = 0, 0.8, 0.9, 0.89, oo for set B. The sum of the Reynolds stresses, magnetic and 
Coriolis forces is shown in black dots and a solid black line with a grey shade showing the 
standard deviation of the time-average. Each panel has a vertical dashed line marking rm- 



6. Discussion 



We find that the geometry and temporal variabihty of the magnetic field 
depends, to first order, on the spatial distribution of force balances. The de- 
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Figure 10: Similar to figure [o] but for set C. From top to bottom increasing values of Xm = 
0, 0.8, 0.9, 0.98, oo. Reynolds stresses, magnetic and Coriolis forces are shown in green, blue 
and red respectively. The sum of them is shown in black dots joined by a solid black line. 



crease of the electrical conductivity at the top of the simulated fluid results in 
significant changes in the resultant flow regime. 



Christensen and Aubert ( 2006 ) reported an increase in the wave number with 



increasing Ra for non-magnetic convective runs, in agreement with experimental 
results (Aubert et al. , 2001). They also point out that for magnetic cases the 
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correlation between mean wave number and Rayleigh number is incoherent, and 
they attribute it to a difference in force balances present in the magnetic cases 
versus the non-magnetic cases. In our models, for weakly forced cases, the 
prograde tilt of the Taylor columns caused by the spherical boundary is allowed 
for the cases where the Lorentz force is small at the boundaries (i.e., Xm = 0-8 
and Xm = 0.9 in set A). Based on the snapshots shown here, e.g.. Figure [6j 
there is indication that the characteristic length scale of the flow depends not 
only on the presence or absence of the Lorentz force but also on its spatial 
distribution. A reduction in the Lorentz force at the top boundary results in 
an overall change in the flow for sets A and C, where the convection becomes 
similar to that of the non-magnetic case. In set B the flow behaviour is different 
for high/low electrical conductivity regions changing the characteristic length 
scale of the flow as a function of radius. 

The models with free-slip boundaries (set C) show that strong azimuthal 
flow in the non-magnetic case becomes weaker in regions of considerable electri- 
cal conductivity (see figure |5|. This result is expected since strong zonal flows, 
can develop for free slip boundaries are damped by magnetic Lorentz forces. 
The Lorentz forces are proportional to the flow velocity and to the strength and 
volume of the conductive layer. In sets other than C, the resultant velocity is 
weaker. This is also expected, since strong zonal flows do not develop in no-slip 
boundary systems. 



6.1. Boundary layers 

It has been suggested in a recent paper that rotating convection can be 



controlled largely by the boundary layers (King et al. 2009). They compare 
the relation between: 1) the thermal boundary layer, 6k ^ where Nu = 

^ p^J^^rj. is the Nusselt number, c is the heat capacity, and Q is the total 

heat flow at the top boundary; and 2) the Ekman boundary layer Se ^ D \fE. 

For non-magnetic rotating convection viscous and Coriolis forces govern the 
dynamics of the Ekman boundary layer. For convective dynamos there also 
exists the magnetic Lorentz force. The non-dimensional magnetic parameter 

analoe^ous to the Ekman number is the Hartmann number B. = . 

which is defined as the ratio of Lorentz to viscous forces. Associated with the 
Hartmann number is the Hartmann boundary layer. Inside the Hartmann layer 
viscous forces are important, while outside the layer magnetic forces dominate. 
The Hartmann layer thickn ess is Sh = y/ji^Xpy^ where is the magnetic 
field normal to the surface (IPotherat et al. 2002 ) . Written in terms of the non- 
dimensional parameters used here, 



where do is a measure of the normalized electrical conductivity in the layer. 

For this study the boundary layers are also affected by radially variable 
electrical conductivity. Decreasing the conductivity near the outer boundary 
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decreases the Lorentz force. A decreased electrical conductivity allows for a 
thicker boundary layer in which turbulent and/or molecular viscous forces are 
important. Figure [TT] shows the time averaged boundary layer thickness of 
Hartmann, Ekman and thermal boundary layers as a function of Xm- In sets 




Figure 11: This figure shows the time-averaged boundary layer thickness in units of as a 
function of Xm- Triangles correspond to 6^ = ^ n (^^^^ time variance is small and is 

not plotted), and circles correspond to 6h = D ^ ^m b- Horizontal dashed lines for the 

Ekman boundary layer thickness are also plotted for sets A, C, and C where 5e = 10~^ and 
for set B where 5e = VlO~^. 

A, C and C the Ekman boundary layer thickness is 5e = 10~^. When Lorentz 
forces do not influence the flow Sh > Se- This is case for Xj^i = 0.8 in sets A, 
C and C as well as and Xm = 0.9 in sets A and C (see also Table [T]). In the 
case of set B where Se ^ 0.0032, Sh > Se only for Xm = 0.8, but for Xm = 0.9 
Sh = 0.001160 is only a factor of 3 times smaller than Se- The transition in the 
flow behavior, found to be a consequence of the force balances at the top of the 
core, is correlated to the relation between the Ekman and Hartmann boundary 
layers. 

When pressure gradients and gravity forces are balanced by the Coriolis force 
the flow is in a geostrophic balance. If the Coriolis is balanced by Lorentz forces, 
the flow is said to be in a magnetostrophic balance. Noting that we keep the 
thermal boundary layer constant for each set, and that S^ is always large (see 
Table [l] and fi gure Ilk, we found that the relative thickness of Hartmann and 



Ekman layers results in force balances at the top of the core allowing for Lorentz 
or Coriolis forces to control the flow, resulting in magnetostrophic {Sh < Se) or 
geostrophic {Se < Sh) systems. 
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Inspection of the radial magnetic field (Figure [8| , reveals the interplay be- 
tween the Ekman and Hartmann boundary layers, and the variable conductivity 
layers near the outer boundary for the various cases. For set A, a relatively low 
Ra results in relatively weak magnetic field. Here the Hartmann layer thick- 
ness is similar to that of the low conductivity layer. For this set, we interpret 
the presence of strong surface radial magnetic field patches near the magnetic 
tangent cylinder to be indicative of the flow fleld, which tends to be deflected 



by a Hartmann layer topography (Potherat et al. , 2002). For Set B, lower E 



and higher Ra^ result in higher A and a 5h that is thinner than the variable 
conductivity layer in the cases with Xm = 0.8 and Xm = 0.9. This results in high 
radial magnetic fleld patches at the outer surface outside the magnetic tangent 
cylinder (i.e., at lower latitudes than the intersection of the tangent cylinder 
with the outer surface) . 

In summary, for thick Hartmann layers, 5h > ^e^ the flow is unaffected by 
the presence of magnetic forces and it may be controlled by the relation between 



Ekman and thermal boundary layers ( King et al.[ |2009|. However for a thin 



layer, 5h ^ ^e^ the Lorentz forces dominate outside the Hartmann boundary 
layer while inside both, thermal and Ekman boundary layers. In our models, 
Ekman and Hartmann boundary layers deflne which force (Coriolis or Lorentz) 
is dominant just inside the region in which viscous forces are signiflcant, so that 
the system follows a geostrophic or a magnetostrophic regime. 



For rotating-convection, iKing et al. (2009) argue that for stress- free bound 



aries, there is a thermal Ekman layer which takes the place of the Ekman layer. 
Such layer would deflne the volume over which the rotation responds to density 
perturbations, modifying plume formation. They argue that the Ekman layer 
is dynamically important in stress-free systems. We flnd a transition in set C 
for A^^^ to be between Xm = 0-9 and Xm = 0.98 This transition coincides 
with the change in the relation between Sh and the traditional Se = ^/E (see 



flgure 11). This is consistent with the relevance of 5e in the dynamics with or 



without stress-free boundary conditions. 

6.2. Planetary implications 

Computational limitations preclude the use of accurate diffusion coefficients 
in modeling the Earth and planets. Furthermore, since we have carried out only 
a few cases that show the correspondence of strong and a weak external magnetic 
flelds to boundary layers with different viscosities and electrical conductivity 
proflles, we are not in a position to attempt to scale our present results to 
conditions in planetary interiors. However, we can come to some preliminary 
conclusions using estimated planetary values of the Ekman and Hartmann layer 
thicknesses. Considering turbulent molecular viscous diffusion we can estimate 
the earth-like Ekman numbers of 10~^ and 10~^^, respectively, yielding Ekman 
layer thickness estimates of order 200 m and 0.2 m, respectively. Considering 



an Elsasser number at the top of Earth's core of 0.3 (Stevenson 2003), the 
Hartmann to Ekman layer thickness of Earth may be estimated to be Sh/Se ^ 
2 where both, Ekman and Hartmann layers thickness are comparable, and a 
magnetostrophic core is expected. 



24 



Observations of Mercury's magnetic field by Mariner 10 and recent MES- 
SENGER flybys are consistent with a Hermean magnetic field that is weaker 
by about two orders of magnitude, but roughly similar in morphology to that 
of Earth (a moderately tilted axial dipole). Because Mercury is small, and its 
rotation rate is quite low (its sidereal day is 58.6 Earth days), its Ekman number 
is about 100 times greater than that of Earth. The Elsasser number scales like 
A (X jVt^ Mercury's A is less than that of Earth by a factor of roughly 1/500. 
Since the Hartmann layer thickness scales like ^ ^e/ a/A we obtain estimates 
for a Hermean Hartmann to Ekman layer thickness to be 5h/Se ^ 100 where 
one expects a non-magnetostrophic flow. 

7. Conclusions 

In our models, the dynamo- generated magnetic field morphology and in- 
tensity are strongly affected by the relative strength of viscous, Coriolis and 
Lorentz forces near the outer boundary. The relative scale of Ekman and Hart- 
mann boundary layer thicknesses is determined by an electrical conductivity 
gradient. For uniform electrical conductivity models Coriolis and Lorentz forces 
are typically of of the same order, yielding a magnetostrophic balance. Low elec- 
trically conductivity near the top boundary can separate the Ekman-Hartman 
layer into a thin Ekman layer and a thicker Hartmann layer, resulting in changes 
of the detailed magnetic field morphology, and a large-scale external magnetic 
field that is relatively weak (see Table [T]). Given the likelihood of a high con- 
centration of light element in Mercury's core, it is plausible that a low electrical 
conductivity layer is present near the core mantle boundary. This means that 
the dynamics of the boundary layers obtained in our models may be applicable 
to conditions in Mercury's core. Our results imply that a low conductivity layer 
is consistent with Mercury's weak observable magnetic field. However, while 
radially variable electrical conductivity is the mechanism studied in this paper, 
it is one of several models that can result in weak magnetic fields. With the 
anticipated arrival of the MESSENGER spacecraft in orbit, we will soon be in a 
position to use detailed mapping to better constrain the relative contributions of 
the various dynamical processes that generate Mercury's global magnetic field. 
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A. Symbols 



Table 2: Symbols used throughout the document 



Symbol 


Description 


units 


a 


polynomial exponent of the electrical conductivity function 




a 


thermal expansion coefficient 




B 


magnetic field induction vector 


(p/ioA,^^)i/2 


c 


heat capacity 


J (kg K)-i 


X 


shell radii ratio 




Xm 


electrically conductive volume radii ratio 




Se 


Ekman boundary layer thickness 


D 


Sh 


Hartmann boundary layer thickness 


D 




thermal boundary layer thickness 


D 


E 


Ekman number 




9 


acceleration of gravity 




go 


acceleration of gravity at the CMB 




H 


Hartmann number 




K 


thermal diffusivity 


m^s~-^ 


^max 


maximum spherical harmonic degree 




A 


model-normalized Elsasser number 




^CMB 


model-normalized CMB Elsasser number 


pfioXift 


^CMB 


model-normalized CMB radial field Elsasser number 


ppoXiCl 


A 


fiuid electrical diffusivity 


2 —1 

m^s ^ 


A 


normalized electrical diffusivity 




A. 


electrical diffusivity at the ICB 


m^s~-^ 


M 


magnetic energy in the fiuid core 






axisymmetric magnetic energy in the fiuid core 




Mo 


magnetic energy at the CMB 




M- 


axisymmetric magnetic energy at the CMB 






magnetic permeability of vacuum 


H m-i 


Nu 


Nusselt number 




V 


kinematic viscosity 


2 —1 

m^s ^ 




Angular momentum 


s-i 


P 


Pressure scalar 


piyQ 




toroidal potential 






poloidal potential 




Prrf 


modified Magnetic Prandtl number 




Pr 


Prandtl number 




Q 


total heat fiux at CMB 






ratio axisymmetric dipole to total magnetic energy 






ratio axisymmetric dipole to total magnetic energy at the CMB 




V 


Unit vector in the radial direction 




p 


Fluid density 


kg m~^ 


Ra 


Rayleigh number 





26 



Table 2: Symbols used throughout the document 



Symbol 


Description 


units 


Rac 


critical Rayleigh number for the onset of convection 




Re 


Reynolds number 




n 


Shell internal radius 


D 




electrical conductivity transition radius 


D 


To 


Shell external radius 


D 


a 


fluid electrical conductivity 


m~^s 


a 


normalized electrical conductivity 


(^i 




electrical conductivity at the ICB 


m~^s 




normalized electrical conductivity at 


(^i 




normalized electrical conductivity at the CMB 




T 


Temperature scalar 


AT 


u 


flow velocity vector 




z 


Unit vector in the direction of the angular momentum 
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